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A multitemporal method to map snow cover in mountainous terrain is proposed to guide Landsat climate 
data record (CDR) development. The Landsat image archive including MSS, TM, and ETM+ imagery was 
used to construct a prototype Landsat snow cover CDR for the interior northwestern United States. Landsat 
snow cover CDRs are designed to capture snow-covered area (SCA) variability at discrete bi-monthly inter- 
vals that correspond to ground-based snow telemetry (SNOTEL) snow-water-equivalent (SWE) measure- 
ments. The June 1 bi-monthly interval was selected for initial CDR development, and was based on peak 
snowmelt timing for this mountainous region. Fifty-four Landsat images from 1975 to 2011 were pre- 
processed that included image registration, top-of-the-atmosphere (TOA) reflectance conversion, cloud 
and shadow masking, and topographic normalization. Snow covered pixels were retrieved using the normal- 
ized difference snow index (NDSI) and unsupervised classification, and pixels having greater (less) than 50% 
snow cover were classified presence (absence). A normalized SCA equation was derived to independently es- 
timate SCA given missing image coverage and cloud-shadow contamination. Relative frequency maps of 
missing pixels were assembled to assess whether systematic biases were embedded within this Landsat CDR. 
Our results suggest that it is possible to confidently estimate historical bi-monthly SCA from partially cloudy 
Landsat images. This multitemporal method is intended to guide Landsat CDR development for freshwater- 
scarce regions of the western US to monitor climate-driven changes in mountain snowpack extent. 

© 2013 Elsevier Inc. All rights reserved. 


1. Introduction 

Continental ice sheets, sea ice, permafrost, and hemispheric-scale 
seasonal snow cover play an important role in regulating the Earth's 
radiation balance, and poleward-equatorial latent heat transport dur- 
ing hemispheric cool-seasons (Barry, 2002). Equally important, 
mountain glaciers and seasonal snow cover (extent) in the form of 
mountain snowpack (depth) feed seasonal streamflow and replenish 
hydrological catchments (Barnett et al., 2005; Winther & Hall, 1999). 
Across the arid western US, mountain ranges serve as seasonal water 
towers that hold and release snowpack freshwater resources through 
successive snow accumulation and melt. Snow-fed streamflow con- 
tributes approximately 50-70% to the total western US annual 
water budget (Cayan, 1996). With an automated temporally discrete 
snow telemetry (SNOTEL) snow-water-equivalent (SWE) measure- 
ment network already in place for hydrological forecasting (Serreze 
et al„ 1999), it is quite clear that spatially-explicit, satellite-derived 
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climate data record (CDR) development can augment western US 
mountain snow-covered area (SCA) monitoring on past, present, and 
future timescales. 

Evidence is mounting for alarming declines in western US moun- 
tain snowpack as well as decreasing Northern Hemisphere spring 
snow cover extent (Brown & Robinson, 2011; Derksen & Brown, 
2012; Hamlet et al., 2005; Mote et al., 2005; Pierce et al., 2008). De- 
clining snowpack trends have been linked to warming springtime 
temperatures that trigger earlier and faster snowmelt (Cayan et al., 
2001; Westerling et al., 2006). Knowles et al. (2006) find in their as- 
sessment of precipitation-snow ratios that snowpack in lower eleva- 
tion zones are melting faster with temperature-driven phase changes 
in water. Therefore, blending spatially derived snow cover CDRs, to- 
pographic models, SNOTEL SWE, and instrumental climate observa- 
tions using time-series analysis techniques could advance efforts to 
examine historical trends in mountain snowpack extent, but has not 
been attempted because of lacking satellite CDR availability. 

Over the past four decades, optical remote sensing has provided a 
critically important data source for observing the Earth’s changing 
ciyosphere. Multispectral imagery acquired by Advanced Very-High 
Resolution Radiometer (AVHRR), Moderate Resolution Imaging 
Spectroradiometer (MODIS), Landsat multispectral scanner system 
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(MSS), Landsat Thematic Mapper (TM), and Landsat Enhanced The- 
matic Mapper Plus (ETM + ) platforms at different pixel resolutions 
enable snow cover retrieval over continental regions. Snow exhibits 
high-moderately high reflectances at visible wavelengths and low re- 
flectances at near-infrared wavelengths (Warren, 1982; Wiscombe & 
Warren, 1980). These spectral characteristics allow for pixel level 
snow retrieval and classification at acceptable thematic and spatial 
accuracy (Dozier, 1984; Dozier & Painter, 2004; Hall & Martinec, 
1985). Specific advantages to monitoring snow cover using satellite 
platforms include daily-weekly image acquisitions, broad-scale spa- 
tial coverage over remote mountain and high latitude regions 
(Dozier, 1989; Dozier & Painter, 2004; Hall & Martinec, 1985; 
Rosenthal & Dozier, 1996), and most importantly, the Landsat mission 
provides a fairly continuous historical image record at recurrent in- 
tervals. This paper uses the entire Landsat image archive to map 
mountain SCA since the 1970s. 

A multitemporal method is proposed to guide Landsat snow cover 
CDR development for freshwater-scarce regions of the western US. 
Historically, the cost of digital imagery hindered multitemporal map- 
ping and use of partially cloudy images. With a now freely available 
image archive at the USGS EROS Data Center (http://eros.usgs.gov/) 
from Landsat's MSS, TM, and ETM + mission, the opportunity to con- 
struct CDRs on climatically relevant timescales is now possible with- 
out the high cost. On the other hand, raw data availability first 
requires that sound methods for Landsat CDR development be 
outlined through a concept description and data demonstration. 
This initial step ensures that Landsat snow cover CDR products are 
spatially and temporally accurate, radiometrically consistent, and in- 
teroperable with MODIS snow products from the outset (Hall et al., 
2002b). 

This paper introduces an operational multitemporal method to 
map SCA in mountainous terrain to support Landsat CDR develop- 
ment. A prototype Landsat snow cover CDR for peak snowmelt is 
demonstrated for an interior northwestern US sub-region, namely 
the central Idaho and southwestern Montana mountains. This conti- 
nental sub-region was selected for its Pacific Ocean influenced clima- 
tology, intermountain hydrological basin significance, and high 
topographic relief within the broader northern Rocky Mountain re- 
gion. This study region is also motivated by the need for more basic 
satellite research on continental snow cover patterns and processes 
in mountainous terrain. The paper is organized as follows. Section 
two outlines a multitemporal method for snow cover mapping in 
mountainous terrain using the Landsat image archive. This section 
describes how Landsat snow cover CDRs are derived using a normal- 
ized SCA equation, and assesses whether the prototype Landsat CDR 
contains systematic biases that arise from missing imagery and/or 
cloud-shadow contamination. The Landsat snow cover CDR is then 
evaluated to determine whether scale parameterization influences 
the probability distribution and time domain variance of SCA. Section 
three presents the results obtained from the multitemporal method 
and prototype SCA time-series. Section four discusses current Landsat 
snow cover CDR development including efforts to monitor 
climate-driven changes in mountain snowpack extent. 

2. Multitemporal method for Landsat CDR development 

2.1. Landsat archival imagery and study region description 

Landsat World Reference System (WRS)-l and WRS-2 path-row 
image acquisitions over the prototype study region depend on the 
sensor and year of acquisition (Table 1). For this June 1 snow cover 
CDR, Landsat MSS, TM, ETM+ SLC-on, and ETM+ SLC-off images 
were collected during the Julian dates 143 and 158 for 1973-2011. 
Scan line coverage (SLC) ‘on’ and ‘off images reflect time periods 
for ‘pre’ and ‘post’ sensor malfunctions on Landsat ETM+. Landsat 
ETM + SLC-off images have stripes that contain missing pixels across 


Table 1 

Archival Landsat imagery used for June 1 snow cover CDR development. 


Landsat 

sensor 

Julian bi-monthly 
interval 

Time period 

Path-rows 

Images 
used (n) 

MSS 1-3 

143-158 

1973-1983 

42-29, 43-29 

5 

MSS 4-5 

143-158 

1982-1995 

39-29, 40-29 

4 

TM 4-5 

143-158 

1984-2011 

39-29, 40-29 

33 

ETM + 7 SLC-on 

143-158 

1999-2003 

39-29, 40-29 

4 

ETM + 7 SLC-off 

143-158 

2003-2011 

39-29, 40-29 

8 


the image. SLC-off pixels have no data and are considered missing 
coverage. All images were pre-processed, mosaicked, and classified 
for snow cover using the multitemporal method outlined below. 

The Landsat snow cover CDR demonstrated is a prototype 
time-series of June 1 SCA for the central Idaho and southwestern 
Montana mountains (Fig. 1 ). This June 1 CDR captures SCA variability 
during peak snowmelt for 1975-201 1 (Figs. 2 and 3). Much of the SCA 
on or near June 1 for this region is confined to mid-high elevations — 
low elevation zones have already melted out. The central Idaho and 
southwestern Montana mountains are situated within a continental 
semi-arid climate zone (Mitchell, 1976) with both local and regional 
hydrological significance. These snow-fed rivers and tributaries form 
the headwaters of the Columbia and Missouri River basins and sup- 
port potable water uses, hydropower generation, fishery migration, 
agricultural irrigation, and recreation. 

2.2. Landsat snow cover CDR design 

Seasonal mountain snowpack across the Idaho-Montana region 
begins to accumulate in mid-October and reaches maximum SWE 
around early-mid April (Cayan, 1996). During the winter season (De- 
cember-March), SCA in mountainous terrain more or less remains 
near maximum coverage at high elevations with lower elevations 
exhibiting daily-weekly variability due to short-term meteorological 
conditions. Once spring arrives, snow cover (snowpack) begins to 
successively melt with increasing solar irradiance, decreasing albedo, 
and warming springtime temperatures. Thus, variability in mountain 
SCA during snow accumulation and snowmelt periods has been, is, 
and will continue to be retrieved by each Landsat overpass. 

Landsat snow cover CDR development is limited by the 16-18 day 
repeat schedule, which roughly corresponds to a bi-monthly interval. 
This temporal structure enables the potential to fully use all available 
Landsat images through compositing, especially in areas where image 
overlap allows 8-9 day coverage. Landsat snow cover CDRs are based 
on SNOTEL SWE bi-monthly sampling intervals centered on days 1 
and 15 of each calendar month (Serreze et al., 1999), where Landsat 
images on days 8-22 are assigned to the mid-month interval, and 
days 23-7 mark the end-beginning month interval (Fig. 2). These dis- 
crete bi-monthly intervals allow for SCA to be mapped at seasonal, 
interannual, and decadal timescales across a defined spatial domain. 
This spatial domain involves a latitude-longitude grid that covers an 
explicit geographic area. This Landsat CDR grid design provides a 
workable geographic scale to manage data volume and overlapping 
images, preserve spectral properties, and facilitate statistical integra- 
tion with gridded ground-based instrumental climate records. 

2.3. Image registration 

Geo-referencing images to a target coordinate system requires ei- 
ther ground control points or ancillary spatial data such as a digital el- 
evation model (DEM) to accurately assign pixels to an exact latitude- 
longitude position (Seidel et al., 1983). Most available Landsat images 
from the EROS Data Center have already been processed to a 
standard-level of geometric and terrain accuracy (http://landsat. 
usgs.gov/Landsat_Processing_Details.php). Geometric errors may 
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Fig. 1. Central Idaho-southwestern Montana snow cover CDR region. The black triangles show SNOTEL SWE measurement sites used for bi-monthly snow cover CDR design. 


exist with older MSS images. This requires careful treatment on an 
image basis when assembling multitemporal sequences. 

For this multitemporal method, unregistered images are geo- 
referenced to a geometrically accurate image (i.e., preferably the 
same sensor with minimal cloud coverage) using an automatic 
tie-point algorithm. Automatic tie-points are selected by correlating 
the unregistered image with a geo-referenced image using a distribu- 
tional grid of stratified points with maximum correlation. An affine 
transformation is then applied to generate a Root Mean Squared 
Error (RSME) statistic. A RSME threshold of less than 0.5 (i.e., half 
the pixel resolution) is the geometric accuracy standard for Landsat 
snow cover CDRs. 


2.4. Radiometric calibration 

Radiometric normalization across multiple sensors and different 
dates of imagery is required to preserve spectral-radiometric proper- 
ties and temporal consistency between multitemporal image se- 
quences (Chander et al., 2009; Chavez, 1996; Song et al., 2001). For 
example, sensor radiometric differences exist between Landsat MSS, 
TM, and ETM+ platforms because spectral bandwidths have been 
placed at different wavelengths along the electromagnetic spectrum. 
This influences retrieval of pixel-level spectral radiance across a con- 
tinuous spectrum. Furthermore, each image was acquired on different 
days where atmospheric conditions and illumination angles in 


67-81 82-97 98-112 113-127128-142143-158159-173174-188 



Bi-Monthly SNOTEL SWE Sampling Intervals 

Fig. 2. Bi-monthly Landsat snow cover CDR design. All available Landsat images from the entire image archive are assigned into discrete bi-monthly Julian sampling intervals. The 
thick black line represents a snowpack accumulation-melt curve based on observed mean bi-monthly SWE measurements for 1976-2010. Mean SWE estimates were derived from 
a regionalization of six SNOTEL sites across the study region (http://www.wcc.nrcs.usda.gov/snow). The gray shading indicates the June 1 Landsat snow cover CDR bi-monthly 
interval. 
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USGS Gauge 06016000 Beaverhead River at Barrets, MT 



USGS Gauge 13118700 Little Lost River BL Wet Creek NR Howe, ID 



Fig. 3. Gauged mean monthly discharge, cubic feet per second (cfs) for Little Lost and 
Beaverhead Rivers in central Idaho and southwestern Montana. Historical mean 
monthly discharge estimates for Little Lost (1958-2011) and Beaverhead (1908- 
2011) gauges were summarized for the period of record. Beaverhead discharge esti- 
mates for Nov. -Feb. were not available for the full period of record. Note: discharge es- 
timates are provisional and subject to revision by the USGS. 

mountainous terrain were variable between images (Meyer et al., 
1993). Given these between sensor and multi-date differences, each 
image must be standardized by calibrating to an absolute radiometric 
scale. 

Landsat imagery downloads from the EROS Data Center are delivered 
in a geotiff format that contains raw unprocessed digital numbers (DN 
values). Radiometric calibration first requires that DN values for each 
spectral band be converted to an at-sensor radiance value using Landsat 
sensor-specific calibration coefficients (Chander & Markham, 2003). 
Calibration coefficients ensure that at-sensor radiance measurements 
acquired during the sensor's mission incorporate corrections for mea- 
surement changes, mechanical malfunctions, and instrument deteriora- 
tion (Chander et al., 2007, 2009). Landsat images used for snow cover 
CDR development have been calibrated using the published calibration 
parameter files (http://landsat.usgs.gov/science_calibration.php). 

The next step in radiometric calibration involves converting the 
at-sensor radiance to a planetary reflectance value. This at-sensor reflec- 
tance conversion can either be implemented for the top-of-the-atmo- 
sphere (TOA) or the surface (Chander et al., 2009; Vermote & 
Kotchenova, 2008). Several atmospheric correction methods have been 
developed, and often include either a partial or full correction that relies 
on ancillary data inputs (Chander et al., 2009; Chavez, 1996; Coppin & 
Bauer, 1994; Song et al., 2001; Vermote & Kotchenova, 2008; Vermote 
et al., 1997). These model inputs include spectral band solar irradiance, 
earth-sun distance, sun elevation, optical thickness, and atmospheric 
aerosols, but even so, it remains unclear whether a full atmospheric cor- 
rection actually produces superior results (Song et al., 2001 ). For Landsat 
snow cover CDRs, a TOA conversion following the methods in Chander et 
al. (2009) was applied using a cosine correction to account for diffuse ir- 
radiance. This TOA conversion was decidedly appropriate at a minimum 
for preserving radiometric consistency across multiple Landsat sensors. 

2.5. Cloud masking 

Cloud cover contamination remains a significant challenge for opti- 
cal remote sensing in planetary regions where climatological patterns 
produce transient, frequent, and persistent cloud cover. Since the 


Landsat image record contains partially cloudy images, a robust 
cloud-masking algorithm to remove cloudy pixels is required. Several 
cloud-masking approaches have been developed for specific 
sensor-platforms (Ackerman et al., 1998; Hagolle et al., 2010; Huang 
et al., 2010; Irish et al., 2006) that use multispectral band criteria to ex- 
amine reflectance change and the thermal band to further separate 
cloudy pixels from non-cloudy pixels. Landsat specific algorithms in- 
clude the automated cloud cover assessment (ACCA) algorithm devel- 
oped by Irish et al. (2006) to estimate cloud cover content on Landsat 
ETM + images. More recently, Huang et al. (2010) developed a Landsat 
cloud algorithm that uses a DEM, the thermal band, and cloud-free con- 
fident forest pixels to identify cloudy pixels. Most studies agree that in- 
cluding the thermal band improves algorithm performance (Hagolle et 
al., 2010; Huang et al., 2010; Irish et al„ 2006). 

Although the ACCA is computationally complex and process inten- 
sively by design, this algorithm is able to separate pixels containing 
snow from pixels containing clouds, and distinguish warm clouds 
from cold clouds using spectral signatures and thermal band cloud pop- 
ulation statistics. The ACCA takes advantage of TM and ETM + bands 2- 
6 by engaging a four-step process that first develops spectral and ther- 
mal band cloud signatures, performs thermal band cloud separation, as- 
signs and aggregates cloudy pixels, and fills cloud holes (Irish et al., 
2006). Lack of middle infrared and thermal spectral bands on MSS limits 
automatic detection of clouds, shadows, and discrimination between 
snow and clouds. Clouds and shadows were subjectively removed by vi- 
sually identifying neighborhoods of dark pixels (i.e„ shadows) that 
were adjacent to brightly illuminated pixels (i.e., clouds). 

2.6. Shadow masking 

Cloud shadows are an inevitable outcome of cloud cover contam- 
ination. The spectral characteristics of cloud-shadow pixels are large- 
ly absorptive with dark low reflectance responses across visible, near 
and middle infrared wavelengths. The aerial extent of cloud-shadow 
contamination is contingent on time of year, time of acquisition, lati- 
tudinal position, topographic relief, cloud area, cloud height, and solar 
illumination. Because cloud shadows share a similar spectral response 
with water and heavily shaded mountainous slopes, simple spectral 
reflectance thresholding is not adequate for shadow detection. 
Hutchison et al. (2009) and Huang et al. (2010) have developed a 
cloud-shadow detection method that includes cloud pixel location, 
cloud pixel height, and solar geometry inputs to project shadow 
pixel location. Once shadow pixel locations have been approximated, 
a moving window is engaged to guide the search for shadow pixels 
exhibiting distinct spectral reflectances (Hutchison et al., 2009). 

Applying the solar geometry cloud-shadow detection method first 
requires a cloud mask. Next, a cloud height estimate is computed for 
each cloudy pixel using both daily surface air temperature (i.e., image 
acquisition date) and Landsat band 6 pixel temperature. Daily mean 
surface air temperature observations used for cloud height estimation 
originate from the NCEP 10-meter reanalysis network (Kalnay et al., 
1996; Masek et al., 2006). For cloudy pixels, normalized temperature 
anomalies are calculated by subtracting band 6 temperatures from 
the surface air temperature. Using the normal lapse rate conversion 
of 6.4 °C decrease for every 1000-meter vertical increase, cloud 
heights for each pixel can be approximated (Huang et al., 2010). 
Below the standard pressure level of 700 hPa, temperature change 
does not always obey the normal lapse rate (Minder et al., 2010). In 
these instances, approximated cloud heights may be substantially 
over or underestimated. Also, actual shadow locations are cast from 
cloud bottoms rather than cloud tops; therefore, thermal band optical 
depth may contribute additional uncertainty to cloud height esti- 
mates. To reduce uncertainty, a range of cloud-shadow projections 
for each pixel is generated using a scale factor (i.e., cloud height mul- 
tiplied by 1.5 or 2 for example) in the denominator of the 
cloud-shadow projection equation (Hutchison et al., 2009; Masek et 
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al., 2006). The final step involves using a moving window with a local 
neighborhood operation to find pixels that show a shadow spectral 
response. A 20 by 20 pixel window (Hutchison et al., 2009) is used 
to search for projected shadow pixels with a spectral reflectance of 
<0.10 for the near infrared (band 4) and <0.10 for the middle infrared 
(band 5) wavelengths. If a pixel satisfies this threshold criterion, it is 
flagged as a shadow (Huang et al., 2010). 

2.7. Topographic normalization 

Land surface radiance represents a combined response to direct 
and diffuse solar irradiance that reflects the sun's illumination angle 
and local topographic angle during time of overpass (Dozier & Frew, 
1990; Lu et al., 2008; Meyer et al., 1993; Smith et al., 1980; Teillet 
et al., 1982). This spectral response is wavelength and land cover de- 
pendent (Meyer et al., 1993), and in mountainous terrain, these ef- 
fects are especially pronounced because of steep slopes and directly 
opposing aspects. As a result, images with mountainous terrain 
need to be normalized for both the sun’s zenith angle and local illumi- 
nation angle. Teillet et al. ( 1 982 ) and others suggest that applying this 
correction is particularly important for multitemporal image analysis, 
and for improving classification accuracy. 

Several quantitative methods for reducing topographic effects 
have been proposed that require empirical or semi-empirical correc- 
tion parameters (Meyer et al., 1993). A key parameter of any correc- 
tion method is whether Lambertian flat surface scattering is assumed. 
Smith et al. (1980) found that for Landsat MSS data, a Lambertian as- 
sumption was generally correct only for slopes below 25° and sun ze- 
nith angles less than 45°. As an alternative, Minnaert's correction is 
generally a preferred method, and requires a derivation of an image 
specific k value that measures the degree to which a surface is 
Lambertian. Landsat CDRs include an image-based Minnaert correc- 
tion to reduce topographic effects. Using a 10-meter USGS DEM and 
sun elevation and azimuth parameters contained within the image 
metadata, illumination layers are generated following the equations 
found in Meyer et al. (1993). The local illumination DEM layer is 
then resampled to 30-meter pixels using a nearest neighbor calcula- 
tion. A stratified random sample from slopes greater than 20° are 
then collected from the local illumination layer and the Landsat illu- 
mination image to derive k. Both populations are log transformed, 
and using linear regression, the slope of the regression line (k 
value) is computed for TM and ETM+ spectral bands 1-5 and 7, 
and MSS spectral bands 4-7. Finally, Minnaert's correction is applied 
back to each spectral band using wavelength specific k values. 

2.8. Image mosaicking 

When mapping regions where multiple Landsat images cover a tar- 
get spatial domain, mosaicking is required to obtain a composite image. 
The main objectives for mosaicking are to maximize cloud-free cover- 
age, piece together images from multiple sensors having different 
image properties, and ensure that original spectral-radiometric proper- 
ties have not been compromised. Landsat CDRs hierarchically select im- 
ages having the lowest cloud cover content to serve as the target image. 
Image mosaicking is then executed using a simple computational over- 
lap function. 

Landsat's 16-18 day image acquisition frequency does limit snow 
cover mapping. In some instances, Landsat's coverage may not be able 
to resolve the difference between resident snow cover and transient 
snowfall (Li et al., 2008). This would result in a tendency to inflate 
SCA, but even so, individual SCA anomaly years could be identified 
using statistical tests and ancillary information from SNOTEL SWE, 
daily MODIS SCA maps, and/or preceding bi-monthly Landsat snow 
cover CDRs. Mosaicking Landsat images from different acquisition 
dates, sometimes even five to seven days, does present a caveat for 
snow cover mapping. 


2.9. Snow cover classification 

Snow is highly reflective at visible wavelengths, and is clearly dis- 
tinguishable from other land surface features (Hall & Martinec, 1985; 
Seidel & Martinec, 2004). Clouds are also highly reflective, but diverge 
with snow at middle infrared wavelengths. While middle infrared 
cloud reflectances remain high, reflectances for snow are absorptive 
(Dozier, 1989; Hall et al., 1995). Spectral bands at visible and middle 
infrared wavelengths allow for spectral discrimination between snow 
and clouds (Dozier & Painter, 2004; Hall et al., 1995). Classifying snow 
cover in multispectral imagery is accomplished using either a binary 
(Seidel & Martinec, 2004) or fractional method (Rosenthal & Dozier, 
1996; Salomonson & Appel, 2004). These two approaches stem from 
the normalized difference snow index (NDSI), a spectral band ratio, 
or subpixel unmixing, a simultaneous linear equation that separates 
vegetation, soils, and lithology endmembers from snow. Both 
methods have been shown to yield acceptable snow cover classifica- 
tion accuracy (Hall & Riggs, 2007; Hall et al., 1995, 1998, 2000, 
2001, 2002a; Rittger et al., 2012; Salomonson & Appel, 2004). The 
NDSI uses TM and ETM + bands 2 and 5, and MODIS bands 4 and 6 
(Hall & Riggs, 2011). Pixels with snow cover greater than 50% have 
been found to have a NDSI value greater than 0.4 (Hall & Riggs, 
2011; Hall et al., 1995). Klein et al. (1998) find that snow beneath 
dense forests can degrade the NDSI response. They suggest including 
a normalized difference vegetation index (NDVI) threshold in con- 
junction with NDSI to improve snow cover detection because the 
spectral response of vegetation with snow has low NDVI values. 

2.9.1. TM and ETM + 

Landsat snow cover CDRs from TM and ETM + use NDSI to retrieve 
and binary classify pixel level snow presence (absence) greater (less) 
than 50% snow cover. This threshold-based snow classification scheme 
for Landsat mimics the MODIS method (Riggs et al„ 2006). Landsat 
pixels are classified as snow if NDSI is greater than 0.4, band 2 is greater 
than 0. 1 1 , and band 4 is greater than 0. 1 0. The band 2 threshold aids in 
the discrimination between snow and water, and the band 4 threshold 
prevents dark targets as being erroneously identified as snow (Riggs et 
al., 2006). At this time, the NDVI threshold (Klein et al., 1998) has been 
excluded because of Landsat's fine spatial resolution, and uncertainty 
whether this response inflates the snow signal during the spring melt 
season. Excluding this criterion from the Landsat NDSI classification 
scheme possibly reduces SCA in dense forest and around the 
snow-covered to snow-free transition zone. 

2.9.2. MSS 

MSS snow-cloud separation is not possible using an automated 
approach. As stated earlier, cloud cover was visually identified using 
shadow detection as a guiding criterion, and then cloud and shadow 
pixels were manually removed. Because clouds and shadows have 
been removed during MSS image processing, highly reflective pixels 
indicate snow cover presence. Using an ISO data-clustering algorithm 
on MSS bands, 30 individual classes are initially identified using spec- 
tral band means and standard deviations. Each class is then 
interpreted to be snow or snow-free. Only cluster classes showing a 
clear homogenous snow cover response are selected for classification. 
Pixels having a mixed snow cover signal are difficult to identify be- 
cause of the mixed spectral responses from snow-free surfaces 
(Dozier & Painter, 2004; Rosenthal & Dozier, 1996). 

2.10. Landsat SCA estimation 

Landsat snow cover CDRs are constructed on bi-monthly time- 
scales to estimate SCA from multispectral measurements over a spa- 
tially explicit domain. Every useable Landsat image from the archive 
is pre-processed, classified for snow cover, and input into a Landsat 
image chronology. An inevitable component of CDR construction is 
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that visible Landsat image coverage will vary over time. For a given 
spatial domain, total land surface area (m 2 ) is calculated using a 
10-m USGS DEM. This area calculation is the total amount of land sur- 
face area that snow can cover at any time including inland water bod- 
ies. This rationale forms the basis for a normalized SCA estimate that 
is independent from image coverage. Landsat SCA percent coverage 
for each annual bi-monthly interval is estimated using the following 
equation: 


Normalized SCA = 


visible SCA^m 2 ^/fland surface area^m 2 j 
- I missing visible coverage (m 2 ) 


+cloud— shadow contamination 



This ratio assumes no change in land surface area, only a change in 
visible coverage and SCA. An estimate of visible coverage also accom- 
panies Landsat snow cover CDRs. Normalized SCA estimates can be 
derived from four possible image combinations: cloud-free, partially 
cloudy, partially cloudy with SLC-off, and partially cloudy with miss- 
ing coverage (Fig. 4). 

2. 31. Landsat CDR homogeneity 

Ground-based, upper air, and satellite CDR homogeneity has and re- 
mains a long-standing data development, production, and distribution 
objective (Peterson et al., 1998). This quality control process is 


particularly important for historical records, and ensures that climate 
data users can access, analyze, and interpret specific datasets without 
concern that ‘non-climatic’ factors such as technological change, instru- 
mentation location change, or measurement protocol have influenced 
the observations in time (Peterson et al., 1998). For historical satellite 
retrievals, these ‘non-climatic’ factors can largely be addressed during 
image pre-processing, but for Landsat snow cover CDRs, additional 
quality assurance steps are required to account for missing image cover- 
age and cloud-shadow contamination when estimating normalized 
SCA. 

If systematic biases are embedded within Landsat snow cover 
CDRs, they will originate from either missing image coverage, clouds, 
shadows, or forest canopy effects (Kane et al., 2008). A primary con- 
cern for Landsat snow cover mapping, especially historical imagery, 
is that the accumulated and combined effects of image acquisition 
frequency, visible coverage, and cloud-shadow contamination will in- 
troduce uncontrollable ‘non-climatic factors’ into the SCA estimation. 
Even though snow cover distribution and scaling across mountainous 
terrain have been shown to be possible (Seidel & Martinec, 2004), 
identifying spatially where potential error could emerge is necessary 
for assessing Landsat CDR quality and homogeneity in time. 

Pixel level relative frequency maps of missing image coverage 
were assembled using no coverage, cloud, and shadow masks derived 
during image processing. MSS and TM/ETM + frequency maps were 
assembled separately because of differing pixel resolution and 
methods for constructing cloud and shadow masks. Relative frequency 
maps for MSS and TM/ETM + visible and non-visible (i.e., combined 


A 


B 





Snow-free 
l~~1 Snow 
I | Cloud 
Shadow 

_ Missing Coverage/SLC-off 


Fig. 4. Possible Landsat image combinations for normalized SCA estimation: a) cloud-free SCA map; b) partially cloudy SCA map; c) partially cloudy SCA map with SLC-off; and 
d) partially cloudy SCA map with missing coverage. 
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missing coverage, clouds, and shadows) pixels were tallied across the 
Landsat CDR. 

2.12. Landsat CDR statistical quality 

Time-series data usually originate from equally spaced continuous 
observations in time (i.e., Landsat images), and are derived based on 
discretely measured values (i.e., SCA) at specific time periods 
(Chatfield, 2004). Time-series theory requires that a non-stationary 
stochastic process with periodic and/or seasonal variations are 
transformed into a stationary process, in this case Landsat SCA, 
where there is no systematic change in the mean or variance over 
time (Chatfield, 2004). Landsat snow cover CDRs are designed to cap- 
ture bi-monthly SCA with the motivation to statistically merge SCA 
with other related climatic time-series such as SWE, temperature, 
and precipitation. Before doing so, it is important to evaluate whether 
the proposed normalized SCA equation produces a scale independent, 
normal SCA probability distribution that is stationary in time. 


To quantify that SCA is scale independent in horizontal and verti- 
cal domains, two separate snow cover CDRs have been derived using 
elevation as a standardized control. The first June 1 SCA CDR has been 
constructed for the full spatial domain, and the second has been 
constructed for the alpine domain above 2500 m. Lilliefor’s test for 
normality is used to determine whether SCA has a normal or 
non-normal statistical distribution with an unspecified mean and var- 
iance (Conover, 1980). Student's t-test for correlated samples, and 
simple correlation is used to test whether SCA mean and variance is 
statistically different, and shares statistically significant covariance if 
constrained by scale parameters. 

3. Results 

June 1 SCA estimates were produced for 30 of the 39 possible years. 
The June 1 SCA CDR for central Idaho and southwestern Montana spans 
1975-2011 with missing years in 1973, 1974, 1978-1982, 1987, 
and 1988 because of either missing coverage or cloud-shadow 
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Fig. 5. Landsat June 1 snow cover CDR: a) June 1 SCA estimates for the full spatial domain; b) full spatial domain standard normal SCA estimates using the 30-year SCA mean and 
standard deviation; c) percent visible image coverage for each full domain SCA estimate; d) June 1 SCA estimates for the alpine domain above 2500 m; e) alpine domain standard 
normal SCA estimates using the 30-year SCA mean and standard deviation; and f) percentage of visible image coverage for each alpine domain SCA estimate. Years 1973, 1974, 
1978-1982, 1987, and 1988 are missing SCA estimates and are marked with line ticks. 
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contamination (Fig. 5). Pixel level relative frequency maps show that if 
SCA estimation is biased, it is more influenced by missing image cover- 
age than cloud-shadow contamination (Fig. 6). Each TM/ETM + 
(Fig. 6a) and MSS (Fig. 6b) frequency map reveals areas across the full 
spatial domain having missing coverage that closely resembles adjacent 
Landsat path-row swath coverage. In the MSS map, scanline striping 
and manually edited cloud-shadow cover polygons are evident. 

Lilliefor's test for normality confirmed that the June 1 SCA 
estimates are non-normal even with a sample size of 30 years 
(Table 2). It is important to note that the June 1 SCA estimates for 
1990 and 2002 have been identified as statistical outliers, and may 
be influenced by transient snowfall at higher elevations during the 
melt season. Student's t-test for correlated samples between sepa- 
rately derived full and alpine domain (i.e„ above 2500 m) June 1 
SCA CDRs indicates a statistical difference between the mean and var- 
iance (Table 3). The correlation coefficient between June 1 SCA CDRs 
is 0.97 (p < 0.001, two-tailed). 

4. Discussion 

The Landsat mission and its moderate spatial resolution, un- 
equaled image record length, and acquisition schedule can supple- 
ment operational SCA mapping on bi-monthly to decadal timescales. 
Bi-monthly snow cover CDR design and development is based on 
using multispectral images from multiple sensors to retrieve SCA dur- 
ing bi-monthly SNOTEL SWE measurement intervals collected across 
the western US. However, this multitemporal snow cover mapping 
method does have inherent limitations that need to be fully acknowl- 
edged from the outset. Primary limiting factors include data gaps in 
historical coverage, image acquisition frequency, and the potential 
for transient snowfall during snow accumulation and snowmelt sea- 
sons. Secondary limitations originate from differences between 
sensor-specific spectral and spatial resolutions, as well as frequent 




Fig. 6. Pixel level relative frequency maps of missing image coverage including cloud 
and shadow contamination: a) TM/ETM + missing coverage for 24 years between 
1985and2011; and b) MSS missing coverage for six years between 1973 and 1 992. 


Table 2 

Lilliefor’s normality test for June 1 SCA estimates. 


Landsat 

CDR 

Record 

length* 

Samples 

(n) 

Lilliefor's Critical-value 

p-value 

Decision 

June 1 SCA- 
Full 

Domain 

1975-2011 

30 

0.176 0.158 

0.017 

Reject 

June 1 SCA- 
Alpine 

1975-2011 

30 

0.160 0.158 

0.044 

Reject 


Decision : The null hypothesis is that SCA estimates come from a normal distribution 
with an unspecified mean and variance. 

Record length*: SCA estimates are not available for all years because of missing imagery 
or cloud-shadow contamination. 


cloud cover contamination over mountainous terrain. For MSS im- 
ages, manual removal of clouds and shadows presents a considerable 
challenge for operational snow cover mapping at the global scale, al- 
though regional mapping efforts may be more feasible. While these 
limitations do present considerable roadblocks when processing 
Landsat images for SCA estimation, taking advantage of additional 
image data from MODIS for example, and/or SNOTEL SWE, may 
offer viable options for filling in temporal data gaps, identifying tran- 
sient snow cover anomalies, and increasing cloud-free coverage. 

Radiometric calibration across the Landsat missions is a critical re- 
quirement for CDR development Feng et al. (2012) have developed 
the Landsat-MODIS Consistency Checking System, and present a frame- 
work for operational Landsat surface reflectance production. In terms of 
this multitemporal method, a full atmospheric correction to surface re- 
flectance is not attempted at this time because of operational and com- 
putational constraints locally. Instead, a TOA correction is applied 
that currently excludes optical thickness, atmospheric aerosols, and 
water vapor corrections common to MODIS products (Vermote & 
Kotchenova, 2008). Moving beyond current constraints into operational 
Landsat CDR production would likely yield a ‘more correct image prod- 
uct’ as the benchmark, but regardless, re-processing of image archives 
can be expected in the future given advancements in algorithm devel- 
opment and computing (Dozier & Frew, 2009). At this time, the above 
TOA method at a minimum provides temporal and radiometric consis- 
tency between Landsat sensors and image inputs into Landsat snow 
cover CDRs. The prototype Landsat snow cover CDR presented should 
be accepted as an initial product, but not final in form. 

Landsat normalized SCA estimation is based on the assumption 
that land surface area over a geographically explicit spatial domain 
constrains and normalizes SCA variance over discrete bi-monthly in- 
tervals. The normalized SCA equation proposed also assumes that re- 
liable SCA estimates can be derived independently from the variation 
in visible image coverage. Cloud-free coverage is not a realistic expec- 
tation for Landsat's image archive, image acquisition frequency, or 
imagery over mountainous terrain. If Landsat's historical image ar- 
chive is to be used comprehensively in innovative ways, accepting 
that target land surface variables, in this case mountain SCA, can be 
derived from partially cloudy images is essential. Currently, Landsat 
snow cover CDRs exclude SCA underneath clouds although spatial 
modeling methods for snow cover extrapolation are available using 


Table 3 

Student's two-sample t-test between June 1 SCA estimates. 



Full domain SCA 

Alpine domain SCA 

Total 

n 

30 

30 

60 

Ex 

2.182 

9.902 

12.084 

Ex 2 

0.293 

5.357 

5.650 

ss 

0.134 

2.088 

3.216 

Mean 

0.072 

0.330 

0.201 

Full SCA mean-Alpine SCA mean 

t 

df 

-0.257 

-6.93* 

29 



Two-tailed significance: p < 0.001*. 
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ancillary data, and have shown acceptable performance for images 
with up to 30% cloud cover obstruction (Seidel & Martinec, 2004). 
Once Landsat snow cover CDRs have been developed for multiple 
bi-monthly intervals, it may be possible to explore whether imple- 
mentation of a SCA extrapolation model is operationally feasible, 
and substantially improves SCA estimation. 

The pixel level ‘salt and pepper’ effect on each frequency map sug- 
gests that over time, cloud-shadow contamination is spatially random 
with greater likelihood of non-randomness directly over mountain 
peaks. Systematic missing image coverage is largely an artifact of 
missing Landsat path-row images on historical timescales. More 
broadly conceived, when missing image coverage is concentrated 
over specific mountain ranges, then SCA estimation may have greater 
uncertainty; however, if missing image coverage is concentrated over 
lower elevation snow-free land, then SCA estimation may by less 
influenced by missing image coverage. At this time, Landsat snow 
cover CDR design and development appears to have few systematic 
biases other than data driven constraints such as missing imagery. 
Landsat CDR homogeneity and quality control should remain a high 
research priority because achieving temporal consistency and spatial 
accuracy is critical for climate-based study in mountainous terrain. 

Landsat snow cover CDRs form the basis for monitoring spatiotem- 
poral variability in mountain SCA, and identifying primary climatic con- 
trols on mountain snowpack extent. The SNOTEL network is able to 
resolve SWE variability at discrete locations, which could be helpful 
for distinguishing between temporal variability in SWE and SCA. One 
primary statistical challenge that accompanies Landsat snow cover 
CDR evaluation is obtaining enough SCA estimates over time to satisfy 
the assumptions of normality and parametric statistical modeling and 
inference. It is well known that instrumental ground-based precipita- 
tion observations have non-normal distributions largely attributed to 
spatial variability in precipitation (Jones & Hulme, 1996). Landsat SCA 
estimates are likely to also be non-normal similar to precipitation. 
This may limit identification of a known SCA mean and standard devia- 
tion with the current record. Achieving normality for Landsat CDRs will 
only come from first, adding more temporal SCA estimates to the prob- 
ability distribution over longer-timescales, and second, correcting for 
transient snowfall in SCA estimation. 

A comparison between full and alpine domain SCA CDRs suggests 
that normalized SCA estimates are stationary, but are statistically 
non-normal. SCA CDRs do show statistically significant covariance in 
time even though scale parameterization influences the mean and 
variance. Increasing the spatial extent of the June 1 snow cover CDR 
to include adjacent regions would likely increase the ability to esti- 
mate SCA more continuously during the 1970s and 1980s. The devel- 
opment of longer Landsat snow cover CDRs on past and future 
timescales will offer greater certainty when determining temporally 
significant SCA trends. 

Multitemporal snow cover mapping in mountainous terrain is 
now possible dating back to 1973 using Landsat. Snow cover CDR de- 
velopment provides a historical starting point for monitoring fine 
scale variability in mountain SCA at bi-monthly timescales and be- 
yond. Missing historical Landsat coverage is inherent to the image re- 
cord, and should be considered a data-driven factor, not a data source 
limitation. This multitemporal method is operationally executable 
and grounded in well-validated algorithm selection and performance. 
Refinements to this multitemporal method should be expected in the 
future given algorithm advancement, greater computing capacity, ad- 
ditional imageiy acquired during the Landsat Data Continuity Mission 
(LDCM), and the potential for multi-sensor integration between 
MODIS and Visible Infrared Imaging Radiometer Suite (VIIRS). 

Using Landsat to monitor historical changes and trends in moun- 
tain snowpack extent through snow cover CDRs is innovative, 
warranted, and can support climate variability and change study, 
hydroclimatological modeling and prediction, and freshwater re- 
source projections under climatic warming. Future work will merge 


Landsat snow cover CDRs with SNOTEL SWE and surface temperature 
and precipitation observations to not only identify the temporal cli- 
matic controls on mountain SCA, but also enable long-term versus 
short-term trend detection over the 20th and 21st centuries. This 
multitemporal method is widely applicable to all mountainous re- 
gions where snow cover accumulates and melts annually, and 
Landsat imagery is available. The proposed method is intended to 
guide Landsat snow cover CDR development across the western US. 
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